function matr = matrix_of_boundary_conditions(nz, nx, CONSTS)

    eta   = CONSTS.eta;
    eps_a = CONSTS.eps_a;

    nxy = func_nxy_from_nz(nz, CONSTS);
    [nxy1 nxy2] = func_nxyk_from_nz(nz, CONSTS);
    ny  = func_ny_from_nxy(nz, nx, CONSTS);
    [J1 J2 J1_1 J2_1] = func_J_k(nxy, nx, CONSTS);
    
    [n1 n2] = func_nk_from_nz(nz, CONSTS);
    
    matr = [-nxy1,              -nxy2,               -1,                   0;  ...
            (1i/eta)*n1.*nxy1,  (1i/eta)*n2.*nxy2,    0,                   -1;     ...
            nxy1.*J1,            nxy2.*J2,            ny./(nxy.^2),        (nz.*nx)./(nxy.^2);   ...
            n1.*1i.*J1_1.*nxy1,  n2.*1i.*J2_1.*nxy2,  (nz.*nx)./(nxy.^2),  -eps_a*(ny./(nxy.^2))];
        
end


